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1.  Introduction 


Micro-mechanics-based  material  models  such  as  the  Tonge-Ramesh  material  model1 
have  the  potential  to  connect  processing  details  like  the  residual  flaw  distribution  to 
simulation  capabilities  that  can  be  used  in  the  design  of  next-generation  armor  sys¬ 
tems,  with  the  ultimate  goal  of  reducing  the  cost  and  lead  time  for  developing  new 
capabilities.  These  goals  require  an  understanding  of  the  performance  of  these  ad¬ 
vanced  material  models  under  conditions  relevant  for  simulating  ballistic  impact 
events.  The  papers  that  initially  presented  the  model  demonstrated  that  the  model 
can  be  used  for  impact  events  where  there  is  moderate  deformation  of  the  materials 
using  a  material  point  method  numerical  scheme.1,2 

The  model  may  be  used  for  both  nonpenetrating  impact  events,  such  as  those  inves¬ 
tigated  by  Tonge  and  Ramesh,  and  impact  events  with  very  large  deformations  such 
as  long-rod  penetration.  Since  nonpenetrating  impact  events  were  addressed  in  the 
initial  articles,  this  work  focuses  on  the  behavior  of  the  model  for  events  with  large 
amounts  of  projectile  penetration. 

This  report  begins  by  summarizing  the  Tonge-Ramesh  material  model  including 
both  the  original  2-surface-distortion/compaction  model  for  granular  flow  and  a 
more  recent  single- surface  granular  flow  model.  Following  this  summary,  the  results 
from  the  Southwest  Research  Institute  (SwRI)  effort  incorporating  the  model  into 
the  EPIC  finite  element  code  are  summarized,  including  their  findings  regarding  the 
limitations  of  the  model.  After  presenting  background  on  the  various  iterations  of 
the  Tonge-Ramesh  model,  2  long-rod  test  problem  geometries  are  discussed.  Simu¬ 
lations  of  these  impact  events  using  a  Lagrangian  mesh  with  SPH  particle  insertion 
in  ALE3D  are  used  to  confirm  that  the  model  behaviors  observed  by  SwRI  are  also 
present  in  an  alternative  numerical  framework.  One  of  the  problem  geometries  is 
also  simulated  using  an  Eulerian  mesh  in  addition  to  the  Lagrangian  simulations  to 
assess  the  ability  of  the  model  to  handle  the  errors  associated  with  material  trans¬ 
port.  While  this  report  only  tests  the  updated  model  in  ALE3D,  the  results  should 
also  apply  if  the  updated  model  is  transitioned  to  EPIC. 
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2.  Background 


The  Tonge-Ramesh  material  model1  is  an  isotropic  micro-mechanics-based  damage 
model  designed  for  simulating  ceramic  materials  subjected  to  impact  loading.  The 
original  form  of  the  model1  described  a  material  failure  process  that  was  initially 
dominated  by  micro-crack  growth  through  a  wing-cracking  mechanism  followed 
by  granular  flow  of  a  comminuted  material.  Three  different  material  models  are 
discussed  in  this  report.  The  models  will  be  referenced  using  the  following  terms: 

•  TR1:  This  is  the  original  model. 

•  TRla:  This  model  uses  a  smooth  single  surface  to  define  the  yield  locus  for 
the  comminuted  material. 

•  TR2:  This  model  is  a  simplified  version  of  the  TRla  model  and  is  described 
in  detail  in  Section  3. 

For  completeness,  the  TR1  and  TRla  models  are  described  in  Sections  2.1  and  2.2. 

2.1  Original  TR  Model 

The  model  was  cast  in  a  finite  deformation  framework  with  a  multiplicative  split 
of  the  deformation  gradient  ( F )  into  an  elastic  portion  (Fe)  and  an  inelastic  por¬ 
tion  (F Gp)  such  that  F  —  FeFG  P.  The  model  specifies  the  Kirchoff  stress  as 
a  combination  of  a  deviatoric  contribution  (Tdev)  and  a  hydrostatic  contribution 
( psJeI )•  The  deviatoric  contribution  is  linearly  proportional  to  the  deviatoric  part 


of  be  =  J~2/3FeF Te\ 


(1) 


The  shear  modulus  ( G )  decreases  with  damage  according  to  an  increased  compli¬ 
ance  model: 


(2) 


The  constants  Zr,  Zn,  and  Zc  represent  the  radial,  normal,  and  coupling  compliance 
of  a  single  penny-shaped  crack  in  an  isotropic  elastic  medium.  They  are  functions 
of  only  the  elastic  moduli.  The  coupling  term  was  added  by  Tonge  and  Ramesh,1 


while  the  other  terms  are  based  on  single-crack  work  by  Grechka  and  Kachanov.3 


A  Mie-Griineisen  equation  of  state  is  used  to  define  the  pressure  ps  associated  with 
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undamaged  material  without  porosity 


Ps(J£,0 )  =  pH{Je) 


£o 

2 


(1  -  Je) 


+  Polm  [ec(Jc)  +  cv(9 


Oo)}. 


(3) 


Here  Je  is  the  thermoelastic  volume  change  ratio,  T0  is  the  Griineisen  coefficient, 
p0  is  the  reference  density,  cv  is  the  specific  heat,  0  is  the  temperature,  ec  is  the  cold 
energy  associated  with  compression  to  Je,  and  pH  is  the  Hugoniot  pressure  given 
by 


Pn(Je ) 


w,Cn2(We) 


(1-S(l-Je))2 

PoCl(l  -  J‘) 


J‘  <  1.0 
otherwise. 


(4) 


Damage  is  incorporated  into  the  bulk  modulus  according  to 


K(D)  =  (K-x  +  D  (Zn  +  4 Zc))  1 .  (5) 


Based  on  the  damaged  bulk  modulus,  the  pressure,  including  the  effects  of  damage, 
the  nonlinear  equation  of  state,  and  temperature,  is 


ps(Je,e,D)  = 


K(D) 

An 


ph(J£ 


1  -  y  (1  -  je) 


+  Poro  [ec(Jc)  +  cv6]  )  .  (6) 


The  damage  variable  D  is  computed  from  the  growth  of  a  distribution  of  microc¬ 
racks 

-Wbins 


D  —  ujk  (sfe  +  Ik)3- 


(7) 


k 

Here  ojp-  is  the  number  density  of  flaws  within  a  family  (or  bin)  of  flaws.  These  flaws 
are  represented  by  an  initial  flaw  size,  Sk,  and  have  grown  an  additional  length,  lk. 
The  flaw  growth  rate  depends  on  the  local  stress  intensity  factor,  Kj,  at  the  crack 
tip  according  to 


ar 


Kj  -  K, 


IC 


-V 


(8) 


A'/  -  0.5A7c  ) 

Once  damage  reaches  a  critical  level,  Dc,  the  granular  flow  phase  of  the  material 
model  is  activated. 


In  the  granular  flow  phase,  there  is  a  distortional  yield  surface  written  in  terms  of 
the  Kirchoff  stress  (r)  defined  by 

-  —  B^j  =  0.  (9) 
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/(t)  =  \J  T  dev  :  Tdev  ~Y  +  A 


tr(r 


The  model  assumes  associated  flow  to  determine  the  direction  of  inelastic  deforma¬ 
tion.  The  rate  of  inelastic  deformation  is  determined  by  assuming  linear  viscoplastic 
behavior  parametrized  by  a  relaxation  time  Top.  The  associated  flow  for  the  distor- 
tional  mechanism  introduces  porosity.  To  account  for  pore  collapse,  an  additional 
yield  surface  that  depends  on  the  pressure  (P  =  —  tr(cr)/3),  distension  (JG p),  and 
total  volume  change  ratio  ( J) 


U(p,jcp,j ) 


_ p _ Eh _ exD  ( _ p? -Pa _ 

Pc-Po  Pc-P0  exp  y  2P0(J«P-1) 
(jGp_i)-(JoGp_i)j2(^) 


(JGP  -  j^)) 


JGP  -  1 


P<P0 

P0  <  P  <  Pc  ■  (10) 

Pc<P 


Numerically  these  2  mechanisms  are  solved  sequentially.  The  distortion  mechanism 
is  computed  first,  introducing  porosity,  followed  by  a  compaction  step  where  the 
porosity  is  reduced.  This  sequential  solution  procedure,  while  reasonable  from  an 
operator- splitting  perspective,  introduces  errors  when  large-strain  increments  are 
introduced  for  loading  steps  near  the  intersection  of  the  2  yield  surfaces.  The  source 
of  the  error  is  shown  in  Fig.  1 .  Starting  from  point  A,  which  resides  on  the  distortion 
yield  surface,  an  increment  in  deformation  results  in  a  trial  stress  located  at  point  B. 
From  point  B,  the  return  algorithm  for  the  distortion  mechanism  returns  the  stress  to 
point  C.  The  pressure  at  point  C  can  be  arbitrarily  large  because  it  does  not  account 
for  pore  collapse.  The  pressure  could  be  so  large  that  the  Mie-Griineisen  equation 
of  state  cannot  be  inverted  to  separate  the  elastic  volumetric  deformation  from  the 
inelastic  volume  deformation.  From  point  C,  the  pore  compaction  return  algorithm 
takes  the  stress  from  point  C  to  point  D.  The  correct  return  location  should  be  point 
E.  When  testing  the  material  model  in  the  EPIC  computational  code,4  a  failure  to 
invert  the  Mie-Griineisen  equation  of  state  caused  a  simulation  to  fail  at  an  impact 
velocity  of  4,500  m/s. 
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Hydrostatic  Stress 


Fig.  1  Schematic  of  error  introduced  by  solving  the  distortion  flow  first  (A-B-C)  followed  by 
the  compaction  (C-D).  The  correct  return  location  is  E. 

2.2  TR  Model  with  a  Smooth  Granular  Flow  Yield  Surface 

To  facilitate  comparisons  with  other  ceramic  failure  models,  such  as  Johnson- 
Holmquist-Beissel  (JHB),5  a  unified  yield  surface  model  was  incorporated  to  rep¬ 
resent  both  the  compaction  and  dilatation  of  the  comminuted  material  using  a  single 
surface  in  terms  of  the  Lode  invariants.  This  variant  of  the  model  is  referred  to  as 
TRla. 

For  a  Cauchy  stress  tensor,  a,  with  a  deviatoric  part,  s  =  cr+pl,  the  Lode  invariants 
are 

r  =  \J s  :  s,  (11) 

*  =  ^tr(o-),  (12) 

J  /  o  \  2/3 

and,  sin  30  =  (  —  )  .  (13) 

2  \  J 2  J 
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The  yield  function6  is  written  as 


f(r,z,9,  -fp,evp)  =  Y{9f\r2  -  Ff{z^p)  \Ff(z,  7P)|  Fc(z,  evp) 

a±  —  03 exp  —  /iqV^3 z'j 


i  -  («*>-*»£  if^< 

(K(e”)-X(ep) 


K, 


otherwise 


M7; p)  =  /*i  +  (A*o  -  A*i)  exp  (-//27 P) 

A"(eP  =  Pi  (po  +  7,  (ln(p3  +p2ep  -  |ln(p3  +p2Cp)|)^ 
K(e;)  =  p4x(ep. 


The  plastic  flow  direction  is 


(14) 

(15) 

(16) 

(17) 

(18) 
(19) 


M  —  a(Ndev  +  {3N  iso). 


(20) 


Here  N dev  and  N lso  are  the  deviatoric  and  isotropic  parts  of  yield  surface  normal, 
j3  is  an  input  parameter  (Jj  <  1),  and  a  is  a  normalization  parameter  to  ensure 
that  M  :  M  =  1.  This  granular- flow  model  uses  the  multistage  return  algorithm.7 
For  all  of  the  cases  investigated  in  this  work,  the  Lode  angle  dependence  is  turned 
off  and  T(0)  =  1.  In  this  form  of  the  model,  the  rate  independent  solution  to  the 
plasticity  problem  is  tracked  in  addition  to  the  rate  dependent  solution  to  provide  a 
more  accurate  integration  of  the  viscous  relaxation  model.8 

2.3  Demands  Placed  on  a  Constitutive  Model  by  Multiple  Host  Codes 

Modeling  and  simulation  is  about  providing  some  insight  into  an  event  based  on 
a  reasonable  set  of  assumptions.  Different  tools  for  solving  the  momentum  and 
energy  balance  equations,  will  make  different  assumptions  and  stress  the  mate¬ 
rial  models  in  different  ways.  This  work  focuses  on  the  Tonge-Ramesh  model  in 
ALE3D,  which  provides  options  to  use  a  fully  Lagrangian  description  of  the  mate¬ 
rial  or  introduce  varying  degrees  of  advection  up  to  fully  Eulerian  calculations.  For 
a  material  model,  a  finite  deformation  Lagrangian  formulation  is  the  least  stressing 
condition,  because  the  material  history  is  accurately  tracked.9  In  advecting  calcula¬ 
tions,  the  relative  motion  between  the  material  and  the  mesh  can  corrupt  the  history 
variables  that  are  carried  with  the  material.  Some  mechanism  must  be  introduced  to 
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deal  with  the  massive  material  distortion  that  will  eventually  distort  a  Lagrangian 
mesh  so  severely  that  the  accuracy  is  lost,  elements  can  invert,  and  a  simulation 
will  crash.  This  work  discusses  the  performance  of  the  Tonge-Ramesh  model  in 
a  long-rod  impact  problem  for  both  Eulerian  calculations  and  Lagrangian  calcula¬ 
tions  with  conversion  of  distorted  elements  to  Smoothed  Particle  Hydrodynamics 
(SPH)  particles. 

2.4  Prior  Evaluations  of  the  Model 

SwRI  implemented  the  Tonge-Ramesh  (TR1)  model  into  the  finite  element  code 
EPIC  and  performed  a  similar  evaluation.4  EPIC  is  a  Lagrangian  finite  element  code 
with  various  erosion  algorithms  to  address  highly  distorted  elements.  In  simulations 
of  various  impact  configurations,  they  concluded  that  the  material  model  (TR1)  was 
relatively  robust  but  slow  compared  to  other  models  like  the  JHB10  model.  Of  all 
the  simulation  runs  that  they  performed,  only  an  impact  at  4,500  m/s  failed  to  run 
to  completion. 

Additional  documentation  of  the  TR1  and  TRla  models  is  available  in  the  user 
manual  for  the  material  model.1 1 

3.  Development  of  a  Faster  More  Robust  Version  of  the  Model 
(TR2) 

This  section  presents  a  second  generation  of  the  Tonge-Ramesh  material  model 
(TR2)  with  some  simplifications  and  enhancements  to  ensure  a  more  robust  solu¬ 
tion  especially  in  the  presence  of  advection  errors.  Even  with  these  improvements, 
significant  advection  of  material  prior  to  failure  will  cause  problems  due  to  the 
smearing  of  the  local  flaw  distributions. 

An  initial  analysis  of  the  TR1  and  TRla  models  indicated  that  there  were  4  areas 
where  a  significant  increase  in  execution  speed  or  robustness  may  be  possible  with¬ 
out  compromising  the  physics  of  the  problem.  They  were  in  the  following  areas: 

•  The  use  of  a  finite  deformation  framework,  in  particular  defining  the  devia- 
toric  part  of  the  elastic  stresses  in  terms  of  b 

•  The  inversion  of  the  Mie-Griineisen  equation  of  state 

•  The  granular  flow  return  algorithms 
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•  Tracking  all  9  components  of  symmetric  rank  2  tensors  in  3  dimensions  when 
only  6  components  are  needed 

The  finite  deformation  framework  that  uses  b  to  define  the  deviatoric  elastic  re¬ 
sponse  is  important  for  polymer  materials,  where  large  elastic  strains  are  possible, 
but  may  be  unnecessary  for  ceramic  materials.  By  assuming  constant  moduli  during 
a  timestep,  we  can  eliminate  the  nonlinear  solve  used  to  invert  the  Mie-Griineisen 
equation  of  state  after  computing  the  stress  satisfying  the  projection  back  to  the 
yield  surface. 

3.1  Elastic  Response 

The  model  uses  an  incremental  formulation  where  the  shear  modulus  is  constant 
and  the  bulk  modulus  is  given  by  the  Mie-Griineisen  equation  of  state.  It  is  shown 
here  for  completeness: 


n 

Us 

dUs_ 

dr] 

Ph 

dPh 

dr] 

P 

U  therm 


P 

Q o 

1  -  s  i  T]  -  s2p2  -  S3r]3 
C0  pi  -  2s2r]  -  3s3?72) 

(1  -  sir]  -  s2r]2  -  s3r ;3)2 
PoUsV 

PoU2  +  ‘IPoUj^r, 

Ph  (1  -  0.5T 0/7)  +  PoT0e 
HP 

(1  -  0.5 Tor])  -  0.5T0Ph 

P  therm  T  Prp. 


(21) 

(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


To  simplify  the  implementation,  the  moduli  are  assumed  constant  over  a  time-step. 
The  stress  and  internal  energy  are  integrated  in  time;  the  temperature  is  not  calcu¬ 
lated.  It  can  always  be  estimated  from  the  total  internal  energy  in  a  post  processing 
step.  For  simulations  using  explicit  time-stepping  schemes,  this  is  a  reasonable  ap¬ 
proximation.  Figures  2  and  3  demonstrate  that  while  the  constant  modulus  assump¬ 
tion  does  introduce  some  error  in  to  the  pressure  and  energy  calculations,  the  solu¬ 
tion  converges  linearly  with  the  size  of  the  strain  increment  and  is  the  error  small  for 
strain  increments  of  0.01.  Using  equation-of-state  parameters  that  are  representative 
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of  boron  carbide  and  a  single  element  driver,  the  pressure  at  a  logarithmic  volume 
strain  of  -0.1  was  29.46  GPa,  when  using  1,000  time  increments.  The  associated 
stored  elastic  energy  per  unit  mass  was  0.5607  MJ/kg.  When  using  a  single  time 
increment,  the  computed  pressure  is  23.22  GPa,  and  the  associated  strain  energy 
per  unit  mass  is  0.921  MJ/kg.  This  represents  an  error  of  22%  in  the  pressure  and 
81%  in  the  strain  energy  per  unit  mass  when  a  single  step  is  used  relative  to  us¬ 
ing  1,000  steps  to  obtain  a  logarithmic  volume  strain  of  0.1.  A  logarithmic  volume 
strain  increment  of  0. 1  with  the  associated  increment  in  pressure  of  nearly  30  GPa  is 
very  large  and  one  should  only  expect  a  constitutive  model  to  produce  a  reasonable 
answer;  one  should  not  expect  an  accurate  pressure  or  strain  energy  measure  with 
such  a  large  increment  in  volume  strain. 


Fig.  2  Error  in  the  pressure  calculation  introduced  due  to  a  consent  modulus  assumption  for 
a  Mie-Griineisen  equation  of  state  for  logarithmic  volume  strain  increments  up  to  0.1  and 
logarithmic  volumetric  strains  of  1 


Approved  for  public  release;  distribution  is  unlimited. 


9 


Fig.  3  Error  in  the  energy  calculation  introduced  due  to  a  consent  modulus  assumption  for 
a  Mie-Griineisen  equation  of  state  for  logarithmic  volume  strain  increments  up  to  0.1  and 
logarithmic  volumetric  strains  of  1 


3.2  Volume  Preserving  Inelastic  Flow 

The  model  allows  for  traditional  volume  preserving  plasticity  with  a  linear-hardening 
law  and  linear  rate  dependence.  The  yield  surface  is  defined  by 


/ja(<T-e p)  =  Ikdevll  -  (Hep  +  To).  (29) 

Here  H  is  the  plastic  hardening  (or  softening)  modulus,  ep  is  the  accumulated  plas¬ 
tic  strain  due  to  this  deformation  mechanism,  and  r0  is  the  initial  shear  strength. 
Optionally,  a  relaxation  time  can  be  provided  for  a  Duvaut-Lions  overstress  vis¬ 
coplasticity  model. 

3.3  Damage  Model 

The  model  uses  an  interacting  micro-crack  model  based  on  a  wing-cracking  mech¬ 
anism  and  handles  the  crack  interactions  using  an  Eshelby  type  self-consistent  ap¬ 
proach.1  The  effective  stress  in  the  neighborhood  of  a  crack  cre  depends  on  both  the 
macroscopic  stress  (Gauss  point  or  material  point  stress,  er)  and  the  damage  level. 
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Micro-crack  growth  is  driven  by  the  stress  intensity  factor  (Kr)  at  the  crack  tips, 
which  depends  on  both  a  wedging  force  Fw  and  the  transverse  stress  normal  to  the 
wing  crack  faces 


K, 


sjn  (/  +  0.27s) 


+  0-22  \fn{l  +  sin  {(j))s). 


The  crack  growth  velocity  is  given  by 


(  Kr-Kjc  \ 

\K  r  +  0.5A/c  / 


(30) 


(31) 


Here  vm  is  the  maximum  crack  growth  velocity.  This  performs  the  same  role  as 
Cr/ac  in  Eq.  8  but  avoids  the  need  to  compute  the  Rayleigh  wave  speed  at  each 
time  step. 


The  crack  growth  velocity  is  a  nonlinear  function  of  the  wing  crack  length.  To 
minimize  errors  associated  with  finite- sized  time  steps  and  promote  smooth  crack 
growth,  we  use  the  following  semi-implicit  scheme  to  update  the  crack  length: 


(n+l  T  ^kl(ln- |-i)  & n  i  ) . 


(32) 


The  inputs  crn,  Dn,  and  sn+ \  are  the  stress  at  the  end  of  the  previous  timestep, 
damage  at  the  previous  timestep,  and  starter  flaw  size  (which  is  a  constant);  these 
are  all  known  at  the  start  of  the  crack  growth  calculation.  This  nonlinear  equation 
for  ln+ 1  is  solved  using  the  TOMS  748  algorithm12  implementation  in  the  Boost 
C++  numerical  library  for  each  flaw  family  (or  bin). 

The  default  behavior  of  the  updated  model  uses  micro-macro  scale  consistency  to 
determine  the  amount  of  inelastic  deformation  due  to  the  microcrack  growth.  In¬ 
elastic  deformation  is  allowed  in  the  flow  direction  associated  with  granular  flow 
(M,  defined  in  Eq.  50)  with  a  magnitude  (Ad)  such  that  the  energy  dissipated  is 
equal  to  the  additional  surface  energy  added  per  unit  volume  due  to  micro  crack 
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growth. 


A 'I' micro  =  7T7,  V  Uk  ((S”«  +  >■?'?  ~  (W  +  &))  (33) 

k 

A 'I'  macro  =  :  AXDM  (34) 

<t  =  <rtr  —  C  :  AA dM.  (35) 


The  stress  at  the  end  of  the  damage  increment  depends  on  the  elastic  moduli  (C), 
trial  stress  (0^),  and  the  increment  in  inelastic  deformation  associated  with  damage 
growth  (A A dM).  The  direction  of  inelastic  deformation  (AT)  is  determined  by 
projecting  the  trial  stress  back  to  the  granular  flow  yield  surface.  This  surface  may 
be  much  smaller  than  the  current  stress  causing  damage  growth.  However,  the  return 
algorithm  for  granular  flow  was  designed  to  handle  large  trial  stresses  relative  to 
the  yield  surface.  Equating  the  microscale  energy  dissipation  (ATmicro)  with  the 
macroscale  dissipation  (ATmacro)  results  in  a  quadratic  equation  for  the  magnitude 
of  the  increment  in  granular  flow  associated  with  damage  growth 


AT micr0  —  <X  tr  •  AA/jAT  —  (C  i  AA/)AT)  :  A\dM 

(<t tr  :  AT)  -  J (<rtr  :  AT)2  -  4  (AT  :  C  :  AT)  AT 


AAn  = 


2  (AT  :  C  :  AT) 


(36) 

(37) 


The  legacy  behavior,  where  the  moduli  are  reduced  with  damage  growth,  is  re¬ 
covered  by  setting  the  parameter  GPRefStress<0.  While  reasonable  for  tensile 
loading  within  the  linear  regime,  the  legacy  behavior  is  questionable  under  large 
compressive  loads.  Particularly,  plate  impact  experiments  conducted  with  a  shock- 
release-reshock  loading  have  suggested  that  the  elastic  wave  speed  does  not  change 
in  boron  carbide  after  the  development  of  damage  under  uniaxial  strain  loading.5 
This  new  energy-based  behavior  introduces  stress  relaxation  that  will  appear  to  re¬ 
duce  the  bulk  modulus  in  tension  and  provide  bulking  in  compression. 


3.4  Granular  Flow 

The  2- surface  granular  flow  model  implemented  in  the  original  version  of  the  Tonge- 
Ramesh  model  is  particularly  inefficient  because  the  intermediate  solution  for  the 
pressure  after  the  linear  Drucker-Prager  solve  can  be  unreasonably  large.  The  sin¬ 
gle  surface  model  avoids  this  issue,  but  the  multi-stage  return  algorithm  used  for 
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the  projection  back  to  the  yield  surface  was  not  suitable  for  arbitrarily  large  strain 
increments.13 

3.4.1  Yield  Surface 

The  yield  function6  is  written  in  terms  of  the  Lode  invariants  of  the  stress  r,  z,  and 
6.  History  dependence  is  addressed  by  defining  2  history  variables  that  are  functions 
of  the  inelastic  strain  rate  due  to  granular  flow  (eGp).  These  separate  the  volumetric 
component  of  the  inelastic  strain  (evp  =  tr(eGP )dr)  from  the  deviatoric  component 
(7 p  =  /■  ||dev[eGP]||dr).  The  yield  function  is  similar  to  the  one  used  in  the  TRla 
model.  The  differences  between  this  form  and  the  TRla  form  are  in  the  Lode  angle 
dependence  (r(9)),  and  the  parameter  a3  was  removed  in  favor  of  placing  z0  inside 
the  exponential.  Additionally  the  changes  in  the  use  of  the  history  variables  make 
the  form  of  the  model  slightly  more  general. 


f(r,z,9, 7P,ep)  =  -r2  —  r2(9)Ff(z,  jp)  \Ff(z,  7P)|  Fc(z,  evp) 

Ff(z,  7P)  =  ai  -  ai  exp  ^a2\/3 (z  -  z0 -  a4V3(z  -  zQ) 


1  -  ^  ^ z\ ,  if  y/3z  < 
Fc(z,  ep)  =  ^  (^)-Y(ep)“ 

1  otherwise 

M7p)  =  Ati  +  (ft o  -  fJ>i)  exp  (— /i27p) 

-o(7 p)  =  4  +  (4  ~  zo )  exp  (~H2lp) 


K 


a  2  — 


“  ifai>0 

ai  1 

)  otherwise 


A"(ep)  =  Pi  (^Po  +  ^  (ln(P3  +P2ep  -  |ln(p3  +p2ep|) 

«(0  =P4^(0- 


(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 


A  valid  yield  surface  requires  the  following: 


•  n  >  a4 


•  ai  >  0,  ai  =  0  is  only  valid  if  a4  >  0 


•  0.  4  >  0 
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•  X  <  z0 


•  0  <  P4  <  1 

•  P3  >  0 

•  pi  >  0 


The  yield  surface  under  triaxial  compression  (a i  <  <r2  =  03)  is  illustrated  in  Fig.  4. 
The  yield  surface  parameters  used  to  generate  this  plot  were  chosen  to  illustrate 
the  different  features  of  the  model  and  are  not  representative  of  actual  material 
parameters. 


Fig.  4  Meridional  profile  of  the  granular  material  yield  surface  showing  the  effect  of  important 
input  parameters 
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The  Load  angle  dependence  follows  a  modified  Reuleaux  profile,  which  enables  an 
analytic  return  in  the  9  direction,14 


b  = 


r2e  -  re  T  1 
2  fe  -  1 


a  =  b  —  re 

IT 

£  = - 9  —  arcsin 

6 


^asin(57r/6  +  9)^ 


r(9) 


a2  +  b2  —  2afecos(^). 


(46) 

(47) 

(48) 

(49) 


Here  fe  is  the  ratio  of  the  triaxial  tensile  strength  to  the  triaxial  compressive  strength. 
This  ratio  must  be  greater  than  0.5  and  no  more  than  1.0.  Nonlinear  Drucker-Prager 
behavior  is  recovered  when  fe  =  1.  Figure  5  shows  the  normalized  octahedral 
profile  as  the  parameter  fe  is  changed.  There  is  an  edge  under  triaxial  extension. 


Fig.  5  Normalized  octahedral  profile  of  the  granular  material  yield  surface  showing  the  effect 
of  the  parameter  re 


The  plastic  flow  direction  is 


M  —  a(N  +  f3N  iso). 


(50) 
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Here  (3  is  a  parameter  between  0  and  1  that  controls  the  degree  of  nonassociativity 
in  the  flow  direction,  N  is  the  normal  to  the  yield  surface  (with  isotropic  component 
7Vlso  and  deviatoric  component  A?" dev),  and  a  is  a  normalization  to  satisfy  |  M  \  = 
1.  For  associative  flow  set  (3  —  1 .  The  edge  and  vertex  cases  are  addressed  by  using 
a  closest  point  projection  in  an  appropriately  mapped  energy  space.  This  definition 
of  the  flow  direction  (M )  requires  an  implicit  solution  of  the  plasticity  return.  The 
algorithm  is  discussed  in  Section  3.7.1. 

3.4.2  Rate  Dependence 

The  original  Tonge-Ramesh  material  model  used  a  Duvaut-Lions  viscoplasticity 
model  with  a  constant  relaxation  time  tqp ■  The  linear  viscosity  model  is  a  reason¬ 
able  first  approximation,  but  a  Bagnold  granular  gas  scaling  model1516  is  able  to 
incorporate  additional  physical  processes.  This  scaling  relationship  arises  from  a 
momentum  balance  when  the  granular  flow  is  behaving  like  a  gas.  When  the  pri¬ 
mary  momentum  transfer  mechanism  between  particles  in  the  granular  flow  is  brief 
particle  collisions,  the  shear  strain  rate  7  scales  with  the  relative  particle  density  p, 
particle  size  d,  particle  coefficient  of  restitution  e,  the  solid  material  density  p0,  and 
the  shear  stress  r 

The  model  uses  this  relationship  as  a  guide  and  assumes  a  relationship  between  the 
inelastic  deformation  rate  (which  is  not  necessarily  volume  preserving)  and  the  total 
overstress  (which  is  not  necessarily  purely  deviatoric) 


e  oc 


/l-p1/3\ 

V  P1/3  ) 


6 


e(l  +  e)p0n 


G 


C  1  :  (<T  ~  (Tqs) 

\/HCr  -  “qsW 


(52) 


The  particle  size  (d)  and  the  relative  packing  density  (p)  are  functions  of  the  history 
variables 


d 


flk>2sk  \  X/3 

(  E 

(53) 

\  Afbins  / 

exp(— e„/3). 

(54) 
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To  facilitate  writing  the  constitutive  model  as  a  collection  of  nondimensional  terms, 
we  define  a  reference  stress  (<rref),  reference  relative  density  pref,  and  reference  par¬ 
ticle  size  o?ref.  From  these  reference  quantities,  the  coefficient  of  restitution  of  the 
particles  during  a  collision  and  the  material  density,  a  reference  strain  rate  can  be 
computed 


rief  ~  G 


(55) 


Using  these  reference  quantities,  the  strain  rate  equation  can  be  written  as  the  prod¬ 
uct  of  dimensionless  groups 


Y / 1  -  PV3  Pref3  A  ( C  1  ■  (<T  -  CTqs) 

ret  \^ref/  \  P1/3  1  -  p)lf )  \d  1  /  sjWcr  -  crqs ||/aref 


(56) 


The  Bagnold  scaling  relationship  was  derived  for  granular  gasses1516  and  breaks 
down  as  the  relative  density  of  the  granular  material  approaches  1  (p  — >  1)  to  ad¬ 
dress  this  case  we  only  compute  the  strain  rate  using  a  maximum  value  of  pref  for 
P- 


It  is  computationally  convenient  to  write  the  strain  rate  sensitivity  as  a  Duviont- 
Lions  viscous  relaxation  model 

eGP  =  — C_1  :  (er  -  crqs) 

tGp 

1  =  ( <£rief\  f  (1  f  (4ef  \  /  7  |o~ref| 

rGP  vm;  \  d )  V  vik-^qS 

If  p  >  pref,  the  fraction  involving  the  relative  packing  density  terms  is  evaluated  to 
1.  The  relative  packing  density  scaling  can  be  disabled  by  setting  pref  <  0  in  the 
input.  Similarly,  the  particle  size  scaling  is  disabled  by  setting  drG  <  0  or  running 
the  calculation  is  run  without  damage  growth. 

The  material  model  has  3  options  for  the  rate  dependence  of  the  granular  flow 
process:  rate  independent  (GPRef StrainRate=0),  Bagnold  granular  gas  scal¬ 
ing15-16  (GPRef  StrainRate>0),  and  a  linear  viscosity  model  (GPRef  StrainRate<0). 


(57) 

(58) 
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The  linear  viscosity  model  is 


1  _  G  |  eref  | 

Top  |  CTref  | 


(59) 


3.5  TR2  Material  Model  Parameters 


The  parameters  for  the  TR2  model  are  summarized  in  Table  1. 
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Table  1  Material  parameters  in  TR2  model 


Keyword 

Symbol 

Units 

Meaning 

MGCO 

Co 

L/T 

Bulk  sound  speed 

MGGammaO 

r0 

Grtineisen  constant 

MGS1 

Si 

Linear  term  in  shock  speed  relation 

MGS2 

S2 

Quadratic  term  in  shock  speed  relation 

MGS3 

S3 

Cubic  term  in  shock  speed  relation 

ShearModulus 

G 

P 

Shear  modulus 

rho_orig 

PO 

M/L3 

Reference  density 

FlawDensity 

V 

1/L3 

Number  density  of  flaws 

MinFlawSize 

^min 

L 

Minimum  flaw  size  in  Pareto  distribution 

MaxFlawSize 

Smax 

L 

Maximum  flaw  size  in  Pareto  distribution 

FlawDist Exponent 

a 

Exponent  in  the  Pareto  distribution 

RandomSeed 

integer 

Seed  for  randomization 

RandomMethod 

integer 

Set  randomization  method  less  than  0  to  disable 

BinBias 

Skews  flaw  binning  for  methods  6  and  7 

FlowStress 

TO 

p 

Initial  J2  flow  stress  (set  less  than  0  to  disable) 

Harden ingModulus 

H 

p 

Hardening  modulus  for  J2  flow 

J2RelaxationTime 

T 

Viscous  timescale  for  J2  flow 

KIc 

Kic 

PL1/2 

Critical  stress  intensity  for  crack  growth 

FlawFriction 

P 

Friction  across  the  crack  faces 

FlawAngle 

<t> 

radians 

Flaw  orientation 

FlawGrowthExponent 

7c 

Exponent  for  crack  growth  law 

FlawGrowthVel 

Vm 

L/T 

Maximum  crack  speed 

CriticalDamage 

Dc 

Damage  to  activate  granular  flow 

GPRef StrainRate 

<7ef 

l/T 

Granular  flow  reference  strain  rate 

GPRef Stress 

^ref 

P 

Granular  flow  reference  stress 

GPRefRhoBar 

Pref 

Reference  packing  density 

GPRef Size 

<4rf 

L 

Reference  particle  size 

GFMSmO 

PO 

Initial  low  pressure  slope 

GFMSml 

Pi 

Final  low  pressure  slope 

GFMSm2 

P2 

Slope  change  rate  with  strain 

GFMSpO 

PO 

Crush  curve  parameter 

GFMSpl 

Pi 

P 

Crush  curve  parameter 

GFMSp2 

P2 

Crush  curve  parameter 

GFMSp3 

P3 

Crush  curve  parameter 

GFMSp4 

Pi 

Shear  pressure  coupling  for  pore  collapse 

GFMSal 

ai 

Extrapolation  of  the  <24  slope  to  zo 

GFMSzOO 

OO 

P 

Initial  granular  hydrostatic  tensile  strength 

GFMSzOl 

z1 

z0 

P 

Final  granular  hydrostatic  tensile  strength 

GFMSa4 

<24 

P 

High  pressure  slope 

GFMSBeta 

P 

nonassociativity  for  granular  flow 

GFMSPsi 

re 

Strength  reduction  ratio  for  triaxial  extension 

In  this  new  model  the  parameters  vm,  eief,  crref,  pref,  and  drc(  are  new.  Many  redundant 
parameters  and  switches  were  removed.  The  form  of  the  yield  surface  softening  is 
different  between  TRla  and  TR2;  however,  for  a  granular  flow  response  that  is 
independent  of  the  yield  surfaces  can  be  made  similar 
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(60) 

(61) 


,TR2 


It  "  =  a2 


TRla  TRla 


CL  q 


+  /i 


TRla 


TR2  _  TRla 

(L ^  —  yL6 


a™2  = 


-.TRla 

h 


if  a™la  ^  0 
otherwise 


a™3  =  a™laexp(a™az™2\/3)  +  //Rlaz™V 3. 


(62) 

(63) 


3.6  TR2  History  Variables 


In  addition  to  removing  redundant  parameters  and  switches,  the  new  model  uses 
fewer  history  variables.  The  history  variables  that  are  specific  to  this  material  model 
are  listed  in  Table  2.  The  subscripted  flaw  family  information  is  repeated  for  the 
number  of  flaw  families  that  are  included  in  the  problem.  The  total  number  of  his¬ 
tory  variables  is  7  plus  3  multiplied  by  the  number  of  flaw  families. 

Table  2 

Model  specific  history  variables 

Name 

Symbol  Units 

Meaning 

BulkModulus 

K 

P 

Bulk  modulus 

ShearModulus 

G 

P 

Shear  modulus 

damage 

D 

Damage  variable 

plasticStrain 

ep 

Plastic  strain  due  to  volume  preserving  plasticity 

plasticEnergy 

P 

Energy  absorbed  by  volume  preserving  plasticity 

epsv 

ev 

cp 

Inelastic  volume  strain  due  to  granular  flow 

gam 

7 P 

Inelastic  deviatoric  strain  due  to  granular  flow 

f lawNumber_k 

OJk 

1/L3 

Number  density  of  flaws  for  bin  k 

starterFlawSize_k 

Sk 

L 

Representative  flaw  size  for  bin  k 

wingLength_k 

Ik 

L 

Wing  crack  length  for  bin  k 

In  addition  to  the  model  specific  history  variables  listed  in  Table  2,  the  material 
model  tracks  the  stored  strain  energy  per  unit  mass  (eeias)  and  the  dissipated  energy 
per  unit  mass  (epias).  These  are  integrated  using  a  first-order  accurate  scheme 
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o 

SO 

< 

&  n-\- 1  • 

(64) 

p 

Act  = 

&  n+l  &n 

(65) 

Ao"  iso 

vu 

<1 

3  K  +  2  G 

(66) 

^^elas 

^"n+1  • 

(67) 

p 

^^plas 

^^tot  ^^elas* 

(68) 

Here  crn  is  the  stress  entering  the  material  subroutine  and  Ae  is  the  increment  in 
strain.  The  stored  strain  energy  per  unit  mass  is  stored  in  the  sse  history  variable, 
and  the  dissipated  energy  per  unit  mass  is  stored  in  the  spd  history  variable. 

3.7  Numerical  Implementation 

This  material  model  is  implemented  as  an  Abaqus  User  Material  and  assumes  that 
the  stress  and  strain  increment  coming  into  the  material  model  interface  accounts  for 
objectivity  requirements  through  the  use  of  an  objective  rate  or  the  an  “  unrotated” 
configuration.  Either  approach  can  be  used  because  there  are  no  tensor  internal  state 
variables  in  the  model. 

Since  rotations  are  not  considered  in  the  constitutive  model  and  there  are  no  inter¬ 
nal  tensor  variables,  all  stress  or  strain  tensors  are  symmetric  second-order  tensors 
allowing  significant  improvements  in  associated  tensor  operations.  In  particular, 
specialized  formulas  are  used  to  determine  the  eigen  values,  eigen  projectors,17  and 
inverse  when  needed. 

The  damage  growth  calculation  uses  the  stress  at  the  beginning  of  the  timestep. 
The  volume  preserving  plasticity  mechanism  is  computed  first  using  a  radial  return 
algorithm.  The  output  of  the  volume  preserving  plasticity  calculation  used  as  the 
trial  stress  for  the  granular  flow  calculation  when  granular  flow  is  active.  The  gran¬ 
ular  flow  algorithm  is  more  involved  due  to  the  pressure  and  possible  Lode  angle 
dependence.  The  numerical  approach  is  described  in  the  following  section. 

3.7.1  Rate  Independent  Projection 

The  modified  Reuleaux  profile14  has  a  hydrostatic  tensile  vertex,  edges  in  triaxial 
compression,  and  the  remaining  surfaces  are  smooth  nonplanar  surfaces.  The  input 
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to  the  return  algorithm  is  a  trial  stress  (crtr),  the  Poisson’s  ratio  (//),  the  parameters 
describing  the  current  location  of  the  yield  surface  (ai,  a  9 ,  a4,  re,  X,  and  k),  and  the 
parameter  determining  the  degree  of  nonassociativity  in  the  plastic  flow  direction 
(fi).  From  the  trial  stress,  the  Lode  coordinates  (ztr,  rtr,  and  9  tr)  and  eigen  projectors 
(Pi,  P 2,  and  P3)  are  calculated.  After  this,  the  return  algorithm  checks  the  hy¬ 
drostatic  tensile  vertex,  then  attempts  returning  to  the  smooth  nonplanar  surface.  If 
the  nonplaner  return  causes  the  stress  to  change  sextants  ( 9  <  —  |),  then  the  return 
should  be  to  the  triaxial  compression  edge. 

The  modified  Reulaux  profile  reduces  the  number  of  unknowns  in  the  nonplanar 
stress  return  to  only  the  hydrostatic  pressure 

<-WP) 

(  Ff(z)  a/2 Fc(z)  if  Fc  >  0  and  Ff>  0 
rc{z)  =  < 

I  0  otherwise 

bT(z)  =  \j ( arc(z ))2  +  rl  -  2arc(z)rti  cos(57t/6  +  9tI) 
b(z)  =  brc(z ) 

L2(z)  =  ( bT(z )  -  b(z)f  +  ((z  -  C-tr)  • 


(69) 

(70) 

(71) 

(72) 

(73) 


The  closest  point  return  in  the  energy  norm  is  zqp  where  zCp  minimizes  L2. 14,18  The 
nonlinear  minimization  problem  is  solved  using  Brent’s  algorithm19  implemented 
in  the  Boost  C++  library.  After  computing  the  radial  coordinate  (rep)  and  Lode 
angle  9C p  are  given  by 


£cp  =  arcsin 


rtrsin(57r/6  +  9t , 
br(zcp) 


tcp  —  \J  a2  +  b2  —  2abcos(£) 
rep  =  fcprc(zcp) 


9Cp  =  arccos 


a?  +  r2 


CP 


b2 


2arCp 


5n 

~6' 


(74) 

(75) 

(76) 

(77) 


If  the  computed  Lode  angle  for  this  return  is  outside  of  the  original  sextant  (6cp  < 
— 7t/6),  then  the  return  will  be  to  the  triaxial  compression  edge  (6*CP  =  — 7t/6).  The 
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minimization  problem  for  this  return  case  is 


Ll  =  (rtr  COS  +  Otr'j  ~  rc(z)^  +  ((z  -  C^tr)" 


(78) 


Once  the  Lode  coordinates  are  computed,  the  return  stress  is  reconstructed  from  the 
principal  stresses  and  the  eigen  projectors  of  the  trial  stress  (Pi,  P2,  and  P3) 


^cp  .  rc P  ( sin(6>CP)  .  , 

=  71 +  71 VW  “  “s(  CP) 

=  ((=  -  rCPsin(eCp)y| 

^CP  rev  ( sin(0CP) 

71 +  71  l  WT~  +  0CP 

<rqs  =  ffjPi  +  a2P2  +  cr3P  3. 


(79) 

(80) 

(81) 

(82) 


This  return  does  not  account  for  the  evolution  of  the  internal  state  variables  dur¬ 
ing  the  timestep.  The  nonhardening  return  algorithm  can  be  used  as  the  basis  for 
a  consistent  hardening  solution.18  At  the  end  of  a  timestep  using  a  consistent  algo¬ 
rithm,  the  stress  lies  on  the  yield  surface  that  has  been  updated  based  on  the  strain 
increment  during  the  timestep 


fi^n+li  evn+l’  Tpn+i)  0- 

To  generate  a  consistent  hardening  solution  we  iteratively  solve 


(83) 


y(v)  =  K  -  vKo  =  °-  (84) 

Here  5e%0  is  the  increment  in  computed  using  the  nonhardening  return  algorithm 
and  e£n.  The  quantity  Se?  is  computed  from  the  nonhardening  return  algorithm  after 
updating  the  history  variables  according  to 


=  evn  +  vKo-  (85) 

For  a  hardening  process,  the  solution  to  Eq.  84  is  bounded  between  0  and  1.  A 
known  bounding  range  for  the  solution  enables  solution  algorithms  that  are  guaran¬ 
teed  to  converge.18  We  do  not  try  to  generate  a  consistent  solution  that  accounts  for 
possible  softening  due  to  evolving  7p. 
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Figure  6  illustrates  the  difference  in  convergence  behavior  when  the  consistent  re¬ 
turn  is  used.  These  simulations  are  for  a  single  element  that  starts  with  an  inelastic 
volume  strain  of  0.6,  and  then  the  porosity  is  compressed  out  of  the  material.  There 
is  minimal  error  in  the  pressure  solution  for  the  consistent  algorithm  even  with  a 
strain  increment  of  0.1.  The  faster  algorithm  that  does  not  account  for  the  evolution 
of  the  yield  surface  during  the  timestep  does  not  seem  to  converge  until  a  strain  in¬ 
crement  of  0.001.  Both  algorithms  converge  to  the  same  pressure- volume  response 
at  a  strain  increment  of  0.001  as  shown  in  Fig.  7. 


Fig.  6  Comparison  of  single  element  uniaxial  strain  compression  simulations  from  an  initially 
distended  state  showing  the  effect  of  strain  increment  size 
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Algorithm  Comparison  for  Ae=0.001 
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Fig.  7  Pressure-volume  response  for  an  initially  distended  material  for  uniaxial  strain  com¬ 
pression  with  a  strain  increment  of  0.001  for  the  consistent  algorithm  and  the  faster  nonhard¬ 
ening  algorithm 


3.7.2  Rate  Dependence 

The  numerical  implementation  of  the  granular  flow  viscoplasticity  model  relies  on 
the  ability  to  project  an  arbitrary  trial  stress  (ertr)  onto  the  quasi-static  yield  surface 
(providing  the  value  for  erqs).  Once  the  projection  onto  the  quasi-static  yield  surface 
is  defined,  the  stress  at  the  end  of  the  timestep  is 


&  n+ 1 


At 

tgp 


'  qs 


1+  Ai 

TGP 


(86) 


This  procedure20  to  update  the  stress  only  requires  tracking  the  current  stress  as  a 
history  variable  and  is  more  robust  with  respect  to  advection  errors  (provided  that 
<rqs  can  be  computed  reliably)  than  the  algorithm  that  explicitly  tracked  crqs  as  an 
extra  set  of  history  variables.  When  using  the  Bagnold-type  granular  flow  rate,  the 
relaxation  time  rGp  is  approximated  using  the  trial  stress  as 


1  =  ( G^reA  f(  1  ~  f  dref\  jf  Kf| 

TCP  V  l^ief I  )  (1  -  plf)  )  V  ^  /  V  V  I  K  -  CTqs 


(87) 
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This  approximation  overestimates  the  quantity  \  \cr— crqs||  by  a  factor  of  y7 1  +  Af/rGp 
to  avoid  solving  a  quintic  equation  numerically.  The  extra  cost  of  numerically  solv¬ 
ing  the  quintic  is  not  justified  because  the  backward  Euler  stress  update  algorithm 
is  only  first-order  accurate.20  Using  the  trial  stress  in  the  calculation  instead  of  the 
stress  from  the  previous  timestep  prevents  division  by  zero.  For  global  timesteps 
that  poorly  resolve  the  granular  flow  timescale  Top,  this  algorithm  will  overestimate 
the  rate  of  stress  relaxation.  However,  for  a  decaying  function,  like  the  overstress  in 
Eq.  86,  the  backward  Euler  update  overestimates  the  updated  function  value.  The 
updated  stress  will  be  outside  the  yield  surface  at  the  end  of  the  timestep  avoiding 
saw  tooth  behavior  caused  by  over  estimating  the  increment  in  plastic  strain. 

3.8  Timing  Results 

The  speed  of  the  improved  model  was  measured  by  simulating  a  cube  of  mate¬ 
rial  subjected  to  simple  shear  deformation.  The  cube  was  meshed  with  1,000  el¬ 
ements  (10  elements  per  side)  and  the  simulation  was  run  for  1,000  timesteps 
with  all  displacement  degrees  of  freedom  specified.  This  test  compared  3  differ¬ 
ent  material  models:  the  original  Tonge-Ramesh  implementation  with  a  2-surface 
description  of  the  dilatation  and  compaction  surfaces,  the  single- surface  implemen¬ 
tation  with  the  multistage  return  algorithm,  and  the  improved-single-surface  model. 

All  implementations  were  compiled  without  optimizations.  The  TR1  model  took 
40.2  //s/zonc/cyclc  to  complete  the  simulation.  The  TRla  model  took  16.8  /is/zone/cycle, 
and  the  improved  TR2  model  took  4.0  s/zone/cycle  when  the  consistent  solve  was 
used  and  2.3  //s/zonc/cyclc  when  it  was  not.  These  results  indicate  that  the  improved 
model  is  5-10  times  faster  than  the  original  implementation.  The  timing  results  are 
summarized  in  Table  3.  For  problems  that  are  dominated  by  the  granular  flow  part 
of  the  calculation,  the  TR2  model  is  significantly  faster. 

Table  3  Summary  of  model  evaluation  times  for  1,000  elements  simulated  for  1,000  timesteps 


Model 

Total  model  evaluation  time 

(s) 

TR1 

165 

TRla 

17.6 

TR2 

4.40 

TR2  with  consistency  correction 

4.70 
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The  basic  single- surface  return  algorithm  in  the  TR2  model  seems  to  be  sufficient 
for  most  cases  and  is  faster  than  using  a  consistent  solution.  The  uniaxial  strain 
compaction  problem  used  to  demonstrate  the  difference  between  the  faster  and  the 
consistent  solution  is  expected  to  be  a  worst  case.  This  model  is  intended  to  solve 
problems  where  the  intact  material  fragments,  then  porosity  is  introduced.  In  these 
types  of  loading  there  will  be  some  re-compression  of  the  material,  but  it  is  not 
expected  to  dominate  the  problem. 

4.  Model  Evaluation  Methods 

4.1  Geometry  from  Prior  Experiments 

There  are  experimental  data  from  2  research  groups  on  penetration  of  confined 
boron  carbide  by  high-density,  long-rod  projectiles.21-22  Based  on  these  prior  ex¬ 
periments,  the  following  3  experimental  geometries  were  identified  to  test  the  per¬ 
formance  of  the  material  model  for  impact  velocities  from  1,000  to  5,000  m/s: 

1)  Boron  carbide  cylinder  19  mm  in  diameter  and  39.6  mm  long  encased  in 
4  mm  of  steel  on  the  top,  sides,  and  bottom.  There  was  a  0.07-mm  interference 
fit  between  the  confinement  and  the  boron  carbide. 

2)  Boron  carbide  cylinder  19  mm  in  diameter  and  39.6  mm  long  encased  in 
1  mm  of  steel  on  the  sides  with  4  mm  on  the  top  and  bottom.  There  was  a 
0.07-mm  interference  fit  between  the  confinement  and  the  boron  carbide. 

3)  Boron  carbide  cylinder  16  mm  in  diameter  and  50  mm  long  encased  in  a 
1.59-mm-thick  Ti6-A14  sleeve  with  a  3. 175 -mm  top  plate  and  12.7-mm  bot¬ 
tom  plate.  The  confining  sleeve  had  a  slip  fit  with  the  target  cylinder. 

Geometries  1  and  2  used  a  2-mm-diameter  tungsten  rod  that  was  150  mm  long. 
Geometry  3  used  a  1 1.43-mm-long  tungsten  rod  that  was  0.762  mm  in  diameter. 

4.2  Simulation  Mesh 

The  Lagrangian  mesh  for  geometry  1  is  shown  in  Fig.  8.  The  mesh  for  geometry 
2  is  nearly  identical,  except  elements  are  removed  from  the  outer  circumference 
of  the  cylinder  to  reduce  the  wall  thickness.  There  are  2  cells  across  the  projec¬ 
tile  radius.  In  these  Lagrangian  simulations,  distorted  elements  were  converted  into 
SPH  particles  according  to  history  variables  that  would  indicate  sufficient  loss  of 
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material  strength.  For  steel  or  tungsten  elements,  the  conversion  occurred  after  the 
equivalent  plastic  strain  in  the  element  exceeded  1 .0.  Within  the  ceramic  material, 
elements  were  converted  once  the  inelastic  shear  strain  (7^)  exceeded  0.1. 


Fig.  8  Slice  view  of  the  initial  Lagrangian  mesh  used  for  the  geometry  1  penetration  simula¬ 
tions 


The  initial  mesh  for  the  Eulerian  simulations  of  geometry  1  is  shown  in  Fig.  9.  The 
projectile  is  better  resolved  with  4  elements  across  the  projectile  radius.  This  mesh 
attempts  to  balance  limited  void  material  with  providing  sufficient  space  for  motion 
of  the  target  after  impact. 
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Fig.  9  Representative  mesh  used  for  the  Eulerian  simulations  of  geometry  1  using  4  elements 
across  the  rod  radius 


The  Lagrangian  mesh  for  geometry  3  is  shown  in  Fig.  10.  There  are  2  cells  across 
the  projectile  radius.  Similar  particle  conversion  rules  were  used  for  this  set  of  sim¬ 
ulations  as  for  geometries  1  and  2.  The  short  rod  and  large  impact  velocity  result 
in  conversion  of  the  projectile  material  into  SPH  particles  very  soon  after  the  im¬ 
pact  event  occurs.  The  primary  purpose  of  this  series  of  3  simulations  is  to  test  the 
material  model  at  the  limits  of  its  region  of  applicability. 
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Tungsten 


Fig.  10  Cut  view  of  the  initial  mesh  used  for  Lagrangian  simulations  of  geometry  3 


4.3  Material  Parameters 

The  baseline  material  parameters  for  boron  carbide  are  taken  from  the  original  ap¬ 
plication  of  this  model  to  boron  carbide.2  This  set  of  material  parameters  does  not 
allow  for  volume  preserving  inelastic  deformation  prior  to  full  damage  development 
{D  >  0.125)  and  does  not  limit  the  strength  of  the  fully  damaged  material  when 
sufficient  pressure  is  applied.  Because  the  model  only  allows  Drucker-Prager  plas¬ 
ticity  with  pore  compaction,  at  very  high  pressures,  the  failed  ceramic  can  sustain 
unreasonably  large  deviatoric  stresses.  These  baseline  material  parameters  are  listed 
in  the  TR1  column  of  Table  4.  To  demonstrate  sensitivity  of  the  results  to  the  dam¬ 
aged  material  model,  a  second  set  of  material  parameters  were  defined  based  on  the 
JHB  parameters  for  damaged  boron  carbide  (set  B).23  These  parameters  are  listed 
in  column  TRla  in  Table  4.  Simulations  using  the  TR2  model  used  parameters  that 
were  equivalent  to  the  TRla  parameters. 
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Table  4  Material  parameters  used  in  TR1  and  TRla  model  evaluation 


Parameter  name 


TR1 


TRla 


Parameter  name 


TR1 


TRla 


useDamage 
useGranularPlasticity 
artificial  Viscosity 
BulkModulus 
rho_orig 

HardeningModulus 

J  2RelaxationTime 

FlawDensity 

MinFlawSize 

FlawDistExponent 

RandomSeed 

BinBias 

FlawFriction 

FlawGrowthExponent 

CriticalDamage 

MicroMechPlaneStrain 

DoFlawInteraction 

JLoc 

GPCohesion 

GPPc 

GPPref 

AbsToll 

Maxlter 

GFMSmO 

GFMSm2 

GFMSpl 

GFMSp3 

GFMSal 

GFMSa3 

GFMSPsi 

MGCO 

MGS1 

MGS3 

MGTheta_0 


1 .000e+00 

1. OOOe+OO 

1 .000e+00 

1. OOOe+OO 

0.000e+0() 

0. OOOe+OO 

2.334e+ll 

2.334e+ll 

2.520e+03 

2.520e+03 

N/A 

0.0 

N/A 

0. OOOe+OO 

2.200e+13 

2.200e+13 

1.000e-06 

1.000e-06 

2.600e+00 

2.600e+00 

3.000e+00 

3. OOOe+OO 

-2.000e+00 

-2. OOOe+OO 

6.000e-01 

6.000e-01 

1 .000e+00 

1. OOOe+OO 

1.250e-01 

1.250e-01 

1  .OOOe+OO 

1. OOOe+OO 

1  .OOOe+OO 

1. OOOe+OO 

2.  OOOe+OO 

2. OOOe+OO 

3.000e+06 

N/A 

1.000e+10 

N/A 

1 .000e+08 

N/A 

N/A 

1.0e-12 

N/A 

20 

N/A 

0.001 

N/A 

0. OOOe+OO 

N/A 

1.000e+10 

N/A 

5.000e-03 

N/A 

288.7e6 

N/A 

288. 7e6 

N/A 

1.0 

9.600e+03 

9.600e+03 

9.140e-01 

9.140e-01 

0. OOOe+OO 

0. OOOe+OO 

2.940e+02 

2.940e+02 

usePlasticity 

useOldStress 

artificial  ViscousHeating 

ShearModulus 

FlowStress 

InitialPlasticStrain 

NumCrackFamilies 

FlawDistType 

MaxFlawSize 

RandomizeFlawDist 

RandomizeMethod 

KIc 

Flaw  Angle 

FlawGrowthAlpha 

MaxDamage 

IncInitDamage 

GPTimeConst 

GPGranularSlope 

GPYieldSurfType 

GPJref 

GPS  urfaceTy  pe 

RelToll 

MaxLevels 

GFMSml 

GFMSpO 

GFMSp2 

GFMSp4 

GFMSa2 

GFMSBeta 

GFMSJ3Type 

MGGammaO 

MGS  2 

MGCv 


0. OOOe+OO 

1. OOOe+OO 

1. OOOe+OO 

1. OOOe+OO 

0. OOOe+OO 

0. OOOe+OO 

1.970e+l  1 

1.970e+l  1 

N/A 

10.206e+9 

N/A 

0. OOOe+OO 

2.500e+01 

2.500e+01 

Pareto  (2) 

Pareto  (2) 

2.500e-05 

2.500e-05 

1. OOOe+OO 

1. OOOe+OO 

7. OOOe+OO 

7. OOOe+OO 

2.500e+06 

2.500e+06 

1.044e+00 

1.044e+00 

5. OOOe+OO 

5. OOOe+OO 

1.260e-01 

1.260e-01 

1. OOOe+OO 

1. OOOe+OO 

7.000e-09 

0. OOOe+OO 

8.000e-01 

N/A 

1. OOOe+OO 

N/A 

2. OOOe+OO 

N/A 

0. OOOe+OO 

1. OOOe+OO 

N/A 

1.0e-08 

N/A 

1 

N/A 

0.001 

N/A 

-1.000e-02 

N/A 

1.429e+00 

N/A 

7.500e-01 

N/A 

1.278e-9 

N/A 

1.0 

N/A 

0 

1.280e+00 

1.280e+00 

0. OOOe+OO 

0. OOOe+OO 

9.620e+02 

9.620e+02 

For  geometry  1  and  2,  the  projectile  and  confinement  were  modeled  using  a  Johnson- 
Cook  plasticity  and  failure  model  with  the  parameters  taken  from  the  experimen¬ 
tal  paper.22  The  metals  in  simulations  of  geometry  3  were  also  modeled  using  a 
Johnson-Cook  plasticity  and  failure  model.  The  tungsten  parameters  were  adapted 
from  the  parameters  for  geometries  1  and  2  by  scaling  the  density.  The  other  metals 
were  modeled  using  library  parameters.24 
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5.  Results 


The  majority  of  the  simulations  with  the  baseline  material  model  ran  to  comple¬ 
tion  when  used  with  element  conversion  to  SPH  particles.  The  baseline  material 
model  ran  in  Eulerian  mode,  but  a  convergence  failure  while  inverting  the  Mie- 
Griineisen  equation  of  state  halted  the  simulation  after  only  1  /js.  The  Eulerian  run 
at  1,480  m/s  using  the  TRla  material  parameter  set  ran  to  57  /xs  before  crashing  due 
to  a  “Q-Stop”  error.  These  errors  are  generally  a  symptom  of  some  other  mesh  tan¬ 
gling  or  mixed  zone  issue.  The  log  file  from  this  run  showed  an  invalid  deformation 
gradient  ( F )  in  the  cell  that  caused  the  fatal  error.  The  TR2  version  of  the  model 
discussed  in  Section  3  ran  the  problem  to  completion.  The  geometry  3  problems  ran 
to  completion  for  all  impact  velocities  up  to  4,500  m/s.  The  results  of  the  simula¬ 
tions  are  presented  in  Table  5  with  the  corresponding  depth  of  penetration  reported 
in  the  experimental  papers.21'22  The  material  locations  at  the  end  of  the  Lagrangian 
simulations  of  geometries  1  and  2  are  shown  in  Fig.  11.  The  final  material  locations 
for  the  Lagrangian  simulations  of  geometry  3  are  shown  in  Fig.  12. 


Table  5  Summary  of  simulation  results 


Geometry 

Velocity 

(m/s) 

Material  model 

Simulation  type 

Sim  DOP 

(mm) 

End  time 

(ms) 

Exp  DOP 

(mm) 

1 

1,427 

TR1 

Lagrange 

0 

35 

6 

1 

1,427 

TRla 

Lagrange 

6.9 

58 

6 

1 

1,480 

TR1 

Lagrange 

0 

60 

20 

1 

1,480 

TRla 

Lagrange 

7.9 

60 

20 

1 

1,480 

TR1 

Euler 

0 

1 

20 

1 

1,480 

TRla 

Euler 

39.2 

57 

20 

1 

1,480 

TR2 

Euler 

35.3 

60 

20 

2 

1,517 

TR1 

Lagrange 

6.2 

60 

35 

2 

2,601 

TR1 

Lagrange 

15.0 

25 

35 

3 

1,500 

TR1 

Lagrange 

0 

60 

11 

3 

3,000 

TR1 

Lagrange 

2.3 

60 

22 

3 

4,500 

TR1 

Lagrange 

4.0 

60 

26 

Approved  for  public  release;  distribution  is  unlimited. 


32 


1,427  m/s 
TR1 


1,427  m/s 
TRla 


1,480  m/s 
TR1 


1,480  m/s 
TRla 


1,517  m/s 
TR1 


2,601  m/s 
TR1 


25  |JS 


Fig.  11  Comparison  of  simulation  results  for  the  Lagrangian  simulations  of  geometries  1  and 
2  at  the  end  of  the  simulation 


For  geometry  1,  simulations  were  run  using  both  the  TR1  and  the  TRla  material  pa¬ 
rameter  sets  for  the  1,427-  and  1,480-m/s  impact  velocities.  Experimentally,  a  dwell 
to  penetration  transition  was  observed  between  these  2  velocities.22  When  using  the 
baseline  material  parameters,  there  is  interface  defeat  at  both  of  these  velocities; 
however,  in  calculations  using  the  TRla  material  parameter  set,  there  is  6.9  mm  of 
penetration  for  the  1,427-m/s  impact  and  7.9  mm  for  the  1,480-m/s  impact.  When 
the  numerical  technique  is  changed  from  a  Lagrangian  mesh  with  conversion  to 
SPH  particles  to  a  pure  Eulerian  formulation,  the  penetration  for  the  1,480-m/s  im¬ 
pact  increases  to  39.2  mm.  One  should  not  observe  such  a  large  difference  in  the 
simulation  results  when  the  numerical  technique  is  changed.  Treatment  of  failure 
and  strongly  history  dependent  materials  within  computational  frameworks  that  can 
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simultaneously  handle  arbitrarily  large  shear  deformations  is  an  ongoing  research 
area.  The  meshes  used  for  these  simulations  are  likely  too  coarse  to  accurately  cap¬ 
ture  the  failure  and  localization  processes  in  either  the  Lagrangian  or  Eulerian  sim¬ 
ulations.  However,  the  goal  of  this  report  was  to  document  that  the  TR1  model  and 
TRla  model  could  be  run  in  ALE3D  and  that  a  modified  model  more  suitable  for 
large  simulations  with  advection  had  been  created.  Understanding  the  reason  for 
the  discrepancy  between  the  Lagrangian  and  Eulerian  approaches  is  a  subject  for  a 
different  report.  Similarly,  the  material  parameters  used  for  this  study  were  chosen 
for  consistency,  not  to  provide  the  best  fit  to  the  experimental  observations. 


1.500  m/s 
3,000  m/s 

4.500  m/s 


Fig.  12  Comparison  of  simulation  results  for  the  geometry  3  long  rod  impact  geometry  60  fis 
after  impact 


Geometry  3  was  only  run  using  the  TR1  model  parameters  in  Lagrangian  mode  with 
conversion  to  SPH  particles.  The  purpose  of  the  simulation  was  to  demonstrate  that 
the  model  would  run  at  impact  velocities  up  to  4,500  m/s.  As  shown  in  Fig.  12, 
the  simulations  ran  to  completion  at  60  fis.  In  the  1,500-m/s  simulation  there  was 
interface  defeat.  The  3,000-  and  4,500-m/s  simulations  showed  some  penetration 
into  the  boron  carbide  but  much  less  than  the  experimentally  observed  penetration. 

The  results  of  the  Eulerian  simulations  of  the  1,480-m/s  impact  using  geometry 
1  are  shown  in  Fig.  13.  To  obtain  the  reported  results  for  the  Eulerian  simula- 
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tions  using  the  TRla  model,  the  error  checking  was  disabled  and  the  rate  depen¬ 
dent  granular  flow  was  disabled  and  replaced  with  rate  independent  granular  flow. 
These  limitations  are  addressed  in  Section  3.  Both  the  TRla  and  the  TR2  model 
without  rate  dependence  in  the  granular  flow  provide  similar  results.  Introducing 
rate  dependence  into  the  granular  flow  decreases  the  amount  of  penetration.  The 
TRla  model  could  not  be  run  in  Eulerian  mode  with  the  rate  dependence  because 
the  results  were  severely  corrupted  by  advection  errors.  The  TR2  model  discussed 
in  Section  3  is  significantly  more  robust  and  does  not  have  the  same  limitations. 
The  Eulerian  simulation  results  show  a  significantly  narrower  penetration  channel 
than  the  Lagrangian  simulations  with  conversion  to  SPH  particles.  Flash  X-ray  im¬ 
ages  taken  during  the  experiment  show  a  penetration  channel  that  is  more  consis¬ 
tent  with  the  Eulerian  simulations  than  the  Lagrangian  simulations,  but  the  reason 
for  the  better  agreement  is  unclear.  Eulerian  simulations  using  both  the  consistent 
and  the  faster  version  of  the  TR2  model  produce  similar  results.  Simulations  us¬ 
ing  a  finer  mesh  were  also  run  for  the  Eulerian  case  and  produced  similar  results. 
A  mesh  convergence  study  using  the  Lagrangian  simulations  with  SPH  conversion 
was  not  performed.  The  SPH  conversion  algorithm  proved  to  be  very  computa¬ 
tionally  expensive,  and  the  convergence  study  was  not  justified  based  on  the  SwRI 
work  demonstrating  that  the  TR1  parameter  fit  may  not  be  reasonable  for  impact 
problems  involving  deep  penetration. 
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TRla  model  with  parameters  fit  to  follow  the  JHBB  model  and  rate  independent  granular  flow 


TR2  model  with  parameters  fit  to  follow  the  JHBB  model  and  rate  independent  granular  flow 


Fig.  13  Comparison  of  simulation  results  for  the  Eulerian  simulations  of  the  geometry  1  long 
rod  impact  problem  using  the  TRla  and  TR2  models  with  material  parameters  fit  to  the  JHB 
set  B  boron  carbide  model  and  an  initial  impact  velocity  of  1,480  m/s  50  /./s  into  the  simulation 


6.  Discussion 

These  simulations  demonstrate  that  the  Tonge-Ramesh  material  model  is  imple¬ 
mented  in  ALE3D  and  can  be  used  to  simulate  complex  penetration  problems  up  to 
4,500  m/s.  As  discussed  by  Holmquist  et  al.,4  the  TR1  model  is  relatively  robust,  but 
the  efficiency  could  be  improved.  The  default  model  parameters  do  not  predict  the 
correct  dwell  to  penetration  transition.  The  strength  of  the  damaged  material  is  too 
high  for  the  baseline  case.  This  is  demonstrated  by  the  simulation  under  predicting 
the  penetration  at  1,480  m/s.  The  simulations  using  the  TRla  model  with  parame¬ 
ters  fit  to  the  JHB  set  B  boron  carbide  model  demonstrate  that  reducing  the  strength 
of  the  damaged  material  increases  the  penetration  and  causes  an  over  prediction 
of  the  penetration  depth.  Additionally,  simulations  using  the  TRla  model  with  the 
JHB  set  B  boron  carbide  model  fit  did  not  predict  any  dwell  at  the  lower  1,427-m/s 
impact  velocity  where  it  was  experimentally  observed.  For  both  the  TR1  and  TRla 
simulations  of  geometry  1  at  1,427  m/s,  the  majority  of  the  ceramic  was  damaged 
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and  eligible  for  granular  flow.  The  difference  between  the  2  results  is  due  to  the 
difference  in  the  strength  of  the  fragmented  material.  This  demonstrates  that  the 
TR1  parameters  produce  a  ceramic  that  has  too  little  resistance  to  damage  but  too 
much  resistance  to  deformation  once  it  has  fragmented.  Independent  experiments 
are  needed  to  calibrate  both  the  intact  material  strength  and  the  damaged  material 
strength.  Work  is  underway  to  develop  a  new  set  of  boron  carbide  material  parame¬ 
ters  that  better  fit  the  response  of  the  Pressure  Assisted  Densification  boron  carbide 
selected  as  the  baseline  material  for  the  Materials  in  Extreme  Dynamic  Environ¬ 
ments  Cooperative  Research  Agreement.  It  may  be  possible  to  calibrate  material 
parameters  based  on  the  dwell  to  penetration  transition,  as  was  done  for  the  JHB 
model,23  but  such  an  approach  undermines  the  usefulness  of  the  confined  ceramic 
impact  experiment  as  a  validation  experiment.  One  should  not  use  the  same  experi¬ 
mental  series  for  both  calibration  and  validation. 

The  Eulerian  simulations  demonstrate  that  the  TRla  and  TR2  models  can  be  used 
in  an  Eulerian  calculation.  The  TR2  model  is  faster  and  more  robust  but  may  not  be 
suitable  for  large  elastic  deformations.  Ceramics  are  generally  not  subject  to  large 
elastic  deformations,  and  without  experimental  data  to  base  a  model  on,  there  is 
little  reason  to  retain  the  added  computational  expense  of  the  finite  elastic  defor¬ 
mation  model  used  in  TR1  and  TRla  for  the  elastic  response  of  the  ceramic.  The 
Eulerian  simulations  showed  better  agreement  with  the  experimentally  measured 
penetration  depth  and  the  shape  of  the  penetration  channel  was  captured  better  by 
the  Eulerian  simulations.  In  general,  one  expects  a  Lagrangian  calculation  to  do  a 
better  job  of  capturing  material  interface  effects  and  damage  propagation.  However, 
the  conversion  to  SPH  particles  could  be  eliminating  the  strength  of  the  material  re¬ 
sulting  in  a  wide-shallow  penetration  instead  of  the  deep-narrow  channel  observed 
in  the  experiments. 
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7.  Conclusion 


This  report  has  presented  representative  ALE3D  simulations  using  the  new  Tonge- 
Ramesh  material  model.1  The  model  performs  reasonably  well  under  impact  load¬ 
ing  conditions  up  to  4,500  m/s.  The  model  can  be  significantly  accelerated  by  mak¬ 
ing  some  simplifying  assumptions  about  the  material  behavior  and  then  optimiz¬ 
ing  the  material  model  to  take  advantage  of  these  assumptions.  The  updated  (TR2) 
model  is  suitable  for  both  Lagrangian  and  Eulerian  simulations  and  seems  to  be 
efficient  enough  for  large  simulation  problems. 

Since  the  predictions  of  a  simulation  are  only  as  good  as  the  material  parameters 
used  in  the  simulation,  a  new  set  of  boron  carbide  material  parameters  should  be 
developed  to  provide  the  best  fit  to  existing  experimental  data.  The  baseline  param¬ 
eters  provide  an  intact  material  that  is  too  weak  (damage  grows  too  easily)  and  a 
failed  material  that  is  too  strong. 
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